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Abstract 

We  report  results  of  the  analysis  of  dynamics  of  fluxons  (topological  solitons,  alias 
kinks)  in  the  system  of  three  long  Josephson  junctions  between  three  bulk 
superconductors  which  form  a  prism.  The  system  is  modeled  by  coupled  sine-Gordon 
equations  for  phases  of  the  junctions.  The  relation  for  the  flux  trapped  inside  the  prism 
reduces  the  system  from  three  equations  to  two.  They  include  dissipative  terms,  and  are 
controlled  by  the  frustration  parameter,  which  measures  the  deviation  of  the  trapped  flux 
from  half  a  quantum.  Analyzing  the  effective  potential  of  the  coupled  equations,  we 
identify  different  species  of  kinks  and  nontopological  "bubble"  modes  in  this  system. 
Solutions  for  compound  kinks  and  for  two  types  of  the  bubbles  are  obtained  in  an 
explicit  analytical  form.  Numerical  tests  demonstrate  that  the  compound  kinks  are 
unstable  against  breakup  into  pairs  of  separating  simple  kinks.  The  bubbles  feature 
metastability,  eventually  breaking  up  into  kink-antikink  pairs.  Kinks  which  connect 
different  potential  minima  and  are,  accordingly,  pulled  by  the  potential  difference,  are 
considered  too.  Using  the  momentum-balance  method,  we  predict  the  velocity  at  which 
these  kinks  should  move  in  the  presence  of  the  dissipation,  Numerical  tests 
demonstrate  that  the  analysis  predicts  the  velocity  in  a  virtually  exact  form.  Inelastic 
collisions  between  the  moving  kinks  are  studied  too. 

The  work  has  been  conducted  in  close  collaboration  with  Dr.  Stanford  Yukon  of 
Division  of  Electromagnetic  Technologies,  USAF  Research  Laboratory,  Hanscom 
AFB,  MA  01731,  USA.  In  particular,  I  was  responsible  for  the  analytical  part  of 
the  work,  while  numerical  calculations  have  been  chiefly  carried  out  by  Dr.  Yukon. 
The  following  text  is  a  preprint  of  a  paper  summarizing  results  produced  by  the 
work  on  this  project.  The  paper,  to  be  co-authored  by  S.  Yukon  and  B.  A. 
Malomed,  will  be  submitted  to  Phys.  Rev.  B. 


PACS  numbers:  74.81. Fa;  74.50.+k;  03.75.Lm;  05.45.Yv 
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1. Introduction 

The  dynamics  of  topological  phase  solitons,  which  represent  quanta  of  the  magnetic 
flux,  i.e.,  fluxons,  in  long  Josephson  junctions  (thin  dielectric  layers  separating  two  bulk 
superconductors)  has  been  a  subject  of  many  experimental  and  theoretical  works 
published  in  the  course  of  more  than  30  last  years  -  see,  e.g.,  original  papers  [1-5]  and 
reviews  [6,7],  A  fundamental  model  of  the  fluxon  dynamics  is  based  on  the  sine-Gordon 
equation  with  additional  terms  which  account  for  various  physical  effects,  such  as 
dissipation  and  driving  fields  [1,5].  Similar  extended  versions  of  the  sine-Gordon 
equations  apply  to  the  description  of  a  number  of  other  physical  media,  such  as  charge- 
density  waves  in  chain  conductors,  dislocations  in  crystals,  etc.  [6,8].  In  the  field  of  the 
Josepshon  physics,  a  great  deal  of  attention  was  attracted  to  the  study  of  diverse 
aspects  of  the  fluxon  dynamics  in  coupled  parallel  junctions  [4].  In  particular,  shuttle 
oscillations  of  many-fluxon  complexes  may  find  an  important  application  to  the 
generation  of  coherent  THz  radiation  [9].  Commonly  adopted  models  for  coupled 
junctions  are  based  on  systems  of  sine-Gordon  equations  coupled  by  terms  which 
account  for  the  magnetic  interactions  between  parallel  layers  [4, 6, 7, 9]. 

Most  works  dealing  with  fluxons  in  coupled  junctions  considered  pairs  of  parallel 
junctions  or  multi-layer  stacks  [4,9].  Another  fundamentally  interesting  configuration  is 
the  one  in  which  three  bulk  superconductors  are  put  together  to  build  a  prism,  which 
gives  rise  to  a  symmetric  triangular  set  of  three  Josephson  junctions  emerging  at 
interfaces  between  the  superconductors  [5,10].  An  essential  peculiarity  of  this  setting  is 
that  the  prismatic  structure  may  trap  magnetic  flux,  which  strongly  affects  the  dynamics 
of  phase  solitons  in  the  triangular  Josephson  set.  The  dynamics  features  a  degeneracy 
when  the  trapped  flux  is  exactly  equal  to  half  a  quantum.  Then,  a  small  deviation  of  the 
trapped  flux  from  this  special  value  plays  the  role  of  frustration  which  lifts  the 
degeneracy. 

It  is  relevant  to  mention  that  symmetric  triangular  sets  of  coupled  nonlinear 
dynamical  elements  occur  in  other  areas,  where  specific  features  of  symmetric  and 
broken-symmetry  states  in  such  setting  have  been  studied  in  some  detail.  These 
systems  include  sets  of  three  effective  potential  wells  for  light  beams  in  photonic  media 
[11],  Bose-Hubbard  systems  with  three  sites  [12],  quantum  phase  transitions  and 
Josephson  oscillations  in  sets  of  three  coupled  potential  wells  [13],  and  nonlocal 
interactions  in  dipolar  Bose-Einstein  condensates  (BECs)  loaded  into  a  set  of  three 
potential  wells  [14],  Symmetric  and  asymmetric  states  of  interacting  solitons  in  three- 
core  systems  were  studied  too,  including  a  triangular  configuration  of  linearly  coupled 
parallel  fiber  Bragg  gratings  [15],  coupled  triplets  of  Gross-Pitaevskii  equations  with 
applications  to  the  description  of  three-component  BEC  [16],  and  dissipative  solitons, 
including  vortices,  in  a  symmetric  system  of  three  linearly-coupled  complex  Ginzburg- 
Landau  equations  [17]. 

Our  objective  in  this  work  is  to  develop  a  systematic  analysis  of  the  phase-soliton 
dynamics  in  the  system  of  three  coupled  long  Josephson  junctions,  using  the  model 
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derived  in  Ref.  [5].  In  Section  2,  we  start  by  introducing  the  model  (due  to  the  constraint 
imposed  by  the  trapped  magnetic  flux,  the  respective  system  of  three  coupled  sine- 
Gordon  equations  is  reduced  to  two)  and  analyzing  various  types  of  static  kinks,  which 
may  exist  in  this  system  as  stationary  solutions  connecting  local  minima  or  saddle 
points  of  the  corresponding  effective  potential  in  the  plane  of  two  phase  variables.  The 
static  kinks  (fluxons)  of  two  essentially  different  types  are  possible  in  the  system, 
"simple"  and  "compound"  ones,  which  correspond,  respectively,  to  the  2 n  steps  in  one 
or  both  independent  phases.  A  solution  for  the  compound  kink  is  found  in  an  explicit 
analytical  form  for  arbitrary  values  of  frustration  parameter.  Further,  direct  simulations 
demonstrate  that  the  compound  fluxons  are  unstable:  the  evolution  leads  to  their 
spontaneous  breakup  into  a  pair  of  separating  simple  fluxons,  while  the  simple  fluxons 
themselves  are  completely  stable  topological  solitons.  The  possibility  of  the  splitting  of 
the  compound  kink  into  the  pair  of  simple  ones  can  be  readily  explained  by  comparison 
of  their  energies. 

The  structure  of  the  model  also  admits  solutions  of  another  type,  in  the  form  of 
"bubbles",  i.e.,  soliton-like  modes  built  on  top  of  a  constant  background,  with  zero 
topological  charge.  Two  analytical  solutions  for  the  bubbles  are  found  in  Section  3.  The 
simulations  demonstrate  that  the  bubbles  are  metastable  modes,  which  persist  for  a 
relatively  long  time,  but  eventually  break  up  into  pairs  of  separating  kinks. 

The  same  system  gives  rise  to  kinks  of  a  different  type,  which  connect  backgrounds 
with  different  energy  densities,  which  correspond  to  deeper  and  shallower  potential 
minima,  alias  a  "true  vacuum"  and  a  "false"  one.  In  the  conservative  system,  the 
potential  tilt  makes  the  existence  of  such  kinks  in  a  steady  state  impossible.  However,  in 
the  presence  of  the  dissipation  terms,  the  balance  of  the  driving  potential  force  and 
dissipation-induced  friction  creates  kinks  (alias  shock  waves,  in  this  context)  traveling  at 
the  corresponding  equilibrium  velocity.  In  Section  4,  we  find  this  velocity  in  an  analytical 
form,  using  the  momentum-balance  analysis  [1 ,5].  The  result  is  found  to  be  in  a 
practically  exact  agreement  with  numerical  findings  (and  the  traveling  kinks  are  found  to 
be  completely  stable).  Continuing  the  work  with  the  moving  kinks,  in  Section  5  we  study 
collisions  between  them,  by  means  of  direct  simulations,  which  demonstrate  that  the 
collisions  are  strongly  inelastic;  in  particular,  some  colliding  pairs  annihilate  into 
persistent  quiescent  breathers.  The  work  is  concluded  by  Section  6. 

2.  Kinks  connecting  true  vacuums 

The  system  that  we  investigate  here  is  formed  by  coupling  three  long  Josephson 
junctions  such  that  each  one  shares  a  common  superconductor  with  its  two  neighbors. 
In  the  usual  notation,  the  generalized  sine-Gordon  equations  for  the  phases  of  the  three 
junctions  (jumps  of  the  phases  of  the  wave  function  of  superconducting  electrons  across 
each  junction),  including  the  Ohmic  losses  with  coefficient  a ,  can  be  written  as  [5,10] 
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2  1 

(^)«-(^L+a(^),  +  J7^sin^+J7^[sin^2+sin(^-^2+^)]  =  0, 

2  1 

(^)«-(^L+«(^),  +  Y^sin^+^-[sin^-sin(^-^+^)]  =  0, 


(1) 


where  it  is  assumed  that  the  ratio  of  the  capacitance  per  unit  length  of  the  third  junction 
to  that  of  the  other  two  junctions  is  7,  and,  simultaneously,  the  same  ratio  for  the 
inductance  per  unit  length  is  I/77  (these  assumptions  make  it  possible  to  account  for  the 
possible  difference  in  characteristics  of  the  three  junctions,  while  keeping  equal  Swihart 
velocities  in  them).  System  (1)  contains  only  two  equations,  because,  in  the  limit  of  zero 
transverse  inductance,  the  phase  of  the  third  junction  is  determined  by  the  constraint 
imposed  by  external  flux  <D  trapped  in  the  system,  ^3=$-^2-^,  where 


(j)  =  2nQ>  /  <E>0  =  n  +  8  is  the  phase  corresponding  to  the  external  flux,  <E>0  is  the  flux 

quantum,  and  parameter  8 ,  which  accounts  for  the  deviation  of  the  trapped  flux  from 
half  the  quantum,  determines  the  system's  frustration. 

The  energy  conserved  by  system  (1 )  in  the  absence  of  the  dissipative  terms  ( a  =  0 ) 
is 


TUJ  |  i  1— 

*= =  [{{[(fit + (*  t +*  (Go,- -  & )  )2 + (At + & ): + 7(wx-(A)j2 

-  cos  </\  -  cos  (j)2  +  cos  (<j\  -(/>2-S)-  C(qo)|  dx, 


(2) 


where  C(oo)  stands  for  the  value  of  the  combination  of  the  three  cosines  at  x  =  ±00  (the 
subtraction  of  this  constant  terms  secures  the  convergence  of  the  integrated  energy 
[10]).  Below,  we  will  also  need  an  expression  for  the  momentum  of  system  (1),  which  is 
also  conserved  if  a  =  0 : 


+00 

p  =  -  j  [ W ),  (ri ),  +  (^2  )x  (<4  l+U  O3  )x  \ 

(3) 

=  -  j  K1  ■ +  n) (A  l  (A  )t  +  0  ■ +  7) O2  )x  (A  l~d (A  )x  {<t>2 1  - 7  (</>2 )x  (A  l  }**• 

—00 

The  potential  part  of  energy  (2)  is 

V (^,^2)  =  -cos^  -cos (j)2  —  cos -<j>2  -(/)),  (4) 

which  is  plotted,  in  the  plane  of  in  Fig.  1  for  frustration  ^  =  0.1^- .  When  8  =  0 ,  all 

local  minima  of  the  potential  have  equal  depths,  and  thus  represent  "true  vacuums"  of 
the  system.  In  this  case,  the  stationary  version  of  Eqs.  (4)  gives  rise  to  kink  solutions 
that  connect  these  minima.  The  kinks  at  8  =  0  feature  the  fractional  flux  content,  that 
appears  in  multiples  of  <1>0  /  3  [5].  When  8  *  0 ,  the  potential  minima  that  are  labeled  {A, 
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B,  ....  G}  in  Fig.  1  are  the  true  (lowest-energy)  vacuum  states,  while  shallower  minima, 

with  a  higher  potential  energy,  that  are  labeled  {a,  b . g},  may  be  called  "false  vacuum 

states".  Symbols  //  in  Fig.  1  mark  saddle  points  of  the  potential. 


(Color  online)  Fig.  1.  Contour  plots  of  potential  function  (4)  with  labeled  minima  and 
saddle  points,  for  S  =  0.\x . 

The  values  of  the  coordinates  at  minima  and  saddle  points  of  potential  F(^,^2), 
are  found  from  equations  dV  I  d<j>x2  =  0,  i.e., 

sin  {</)* )  -  sin  {</\  -  ft2  -  c> )  =  0, 
sin  ^(f2 )  +  sin  (j\  -  </>2  -  c) )  =  0, 


which  yields 
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(j\  =(l/3)[(;r  +  £)  +  2;r(«  +  2m)], 

’  (6) 

(j)2  =(1/3)[-(,t  +  <J)  +  2^-(«  +  1)], 

where  n  and  mare  arbitrary  integers.  In  particular,  for  d  =  7rl\0  Eqs.  (6)  give  rise  to 
the  grid  of  minimum  points  displayed  in  Fig.  1 ,  namely,  point  A  for  m  =  1 ,  n  =  0: 

)=-{-n  +  8,n-5)^>  (-0.942, 0.942) ,  (7) 

and  point  b  for  m  =  n  =  0: 

{fi,£)  =  ^(7t  +  S-7c-8)^>(- 1.151,1.151).  (8) 

All  other  minima  may  be  found  by  adding  or  subtracting  multiples  of  In  to  points  (7) 
and  (8).  Further,  the  saddle  point  uAb  marked  in  Fig.  1  has  coordinates  ($*,^*)  =  (£,-£), 

from  which  all  other  saddle  points  can  be  obtained  by  adding  or  subtracting  multiples  of 

In . 


Static  kink  solutions  connect  adjacent  minima  of  the  potential  with  equal  depths, 
shifted  by  In  in  either  or  both  directions  (for  instance,  simple  kinks  may  connect  points 
A  and  C  or  C  and  B  in  Fig.  1,  while  a  "compound"  kink  links  points  A  and  B  in  the 
diagonal  direction).  For  these  solutions,  the  stationary  version  of  Eqs.  (1)  reduces  to  the 
equations  of  motion  for  a  test  particle  with  two  degrees  of  freedom,  ($,^2)  and  unit 

mass,  in  the  pseudo-potential  t7($,^2)  =  -F($,^2)  >  with  coordinate  x  playing  the  role  of 
time.  In  the  rest  of  the  paper,  we  focus  on  the  most  fundamental  case  of  7  =  1  in  Eqs. 
(1). 


For  the  "compound"  kink  solution  connecting  points  A  and  B,  and  passing  through 
saddle  point  uAb  and  the  shallower  minimum  b,  one  can  set  <t\=<j>2.  After  simple 

manipulations  with  the  resulting  mechanical  equation  for  the  single  remaining  degree  of 
freedom,  the  kink  solution,  ^(x) ,  can  be  found  in  an  explicit  form: 


^AB  (x,c>)  =  2tan  1 


AacyxlS-Aa(\c\2 

-arc 

,)e“,W5  -ay 

[2a  | 

[ar  -cr  +a)eaxl2^  -  aj 

8 

<N 

l _ l 

1 - 1 

8 

1 

8 

1 

o' 

1 

V. 

^3 

(9) 


The  constants  here  are  defined  via  roots  ti  2  3A  of  the  quartic  polynomial  which  emerges 
in  the  course  of  the  integration  of  the  equation  for  the  kink's  shape, 


Uq  +  (7| /  +  u2t  +  Cl^t^  +  —  0  , 


(10) 


Nonlinear  Dynamics  of  Globally  Coupled  Sine  Gordon  Equations 


with  t  =  tan(^  /  2) ,  and  the  coefficients  determined  by  values  of  frustration  8  and 
effective  energy,  E  =  6U(</>t)  =  2  cos  fa  -  cos  (2$  -8),  taken  at  the  points  connected  by  the 
kink  solution: 


a0  =  E  -  2  +  cos  8, 
ax  =  4  cos  8, 

a2=2E-6cos8,  (11) 

a3  =  -4  sin  8, 
a4=E +  2  +  cos  8. 

Namely,  ar=Re(tl),t3=cr+ici,  and  a  =  yj(a*  -2arcr+\c\2>Ja4  .  An  example  of  the  kink 

solution  is  displayed  in  Fig.  2  for  5  =  0.1;r  (it  is  relevant  to  mention  that  this  solution  is 
different  from  those  derived  in  Ref.  [18]  for  an  "asymmetric"  double  sine-Gordon 
equation). 


Fig.  2.  (Color  online).  The  plot  of  the  "compound"  kink  connecting  potential  minima  A 
and  B  in  Fig.  1 ,  for  8  =  O.br .  The  red  and  blue  lines  depict  the  completely  coinciding 
analytical  solution  (9)  and  its  numerically  found  counterpart. 

The  next  step  is  to  simulate  the  evolution  of  the  compound  kink,  within  the 
framework  of  full  equations  (1),  to  test  its  stability.  As  shown  in  Figs.  3  and  4,  the  result 
is  that  the  compound  kink  is  unstable  against  numerical  noise,  eventually  splitting  into  a 
pair  of  "simple"  ones,  which  link  pairs  of  points  AB  and  CB  in  Fig.  1  (recall  the 
compound  kink  directly  connects  A  and  B).  Note  that,  unlike  the  exact  solution  (9)  for 
the  compound  kink,  explicit  analytical  solutions  for  the  simple  kinks  are  not  available. 

The  secondary  kinks  produced  by  the  breakup  of  the  compound  one  separate, 
eventually  coming  to  a  halt  due  to  the  action  of  the  friction  terms,  accounted  for  by 
coefficient  a  in  Eqs.  (1).  The  breakup  AB->  AC  +  CB  is  additionally  illustrated  by  Fig.  5, 
which  shows  the  parametric  dependence  between  components  fa  and  in  the  initial 
and  final  states  of  the  phase  fields. 
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(a) 


(b) 
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Fig.  3.  (Color  online)  The  plots  of  </\(x,t)  (a)  and  </>2(x,t)  (b)  for  the  unstable  compound 

kink  of  type  AB,  splitting  into  simple  ones  of  types  AC  and  CB,  for  5  =  OAk  and 
dissipation  coefficient  a  =  0.04 . 


Fig.  5. (Color  online).  Parametric  plots  in  the  plane  of  ($,$, )  illustrating  the  breakup 

of  the  compound  kink,  AB-»  AC  +  CB.  The  blue  straight  line  corresponds  to  the 
stationary  configuration  at  t  =  0;  the  broken  red  line  shows  the  same,  but  at  t  =  200. 

The  possibility  fo  the  breakup  may  be  readily  explained  by  the  consideration  of  the 
system's  energy  (2).  The  substitution  of  the  analytical  solution  (9)  for  the  compound  AB 
kink  solution,  and  of  numerical  solutions  for  the  simple  AC  and  CB  ones,  into  the 
integral  in  Eq.  (2)  yields  the  respective  values  of  the  energy  (rest  masses): 
m ^  =  25.549,  mACCB  =9.527  .  The  difference  between  them  is  a  straightforward 

explanation  to  the  feasibility  of  the  splitting. 

3.  Bubble  modes 

To  get  the  bubble  solutions,  which  seem  as  solitons  set  on  top  of  a  constant 
background,  featuring  identical  values  of  the  phases  at  x  =  ±co  [19],  we  start  from  the 
false  vacuum  state  (shallower  minimum)  b  in  Fig.  1,  and  allow  the  test  particle  to  travel 
in  the  anti-diagonal  direction,  ^  =  -(f>2 .  The  corresponding  potential-energy  profile  is 
shown  in  Fig.  6. 
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Fig.  6.  (Color  online)  The  effective  potential  energy,  U(</>{),  corresponding  to  the  test 
particle  passing  through  point  b  in  Fig.  1  in  the  direction  of  <j>x  =  ~(f>2 ,  for  8  =  0.1;r . 


In  terms  of  Fig.  6,  a  "larger"  bubble  corresponds  to  the  motion  from  point  b 
($  =1.15 )  to  <t>x  =4.81  and  back,  while  the  trajectories  generating  another,  "smaller", 

bubble  bounces  back  from  <j>x  =  -0.21 .  Obviously,  the  two  bubble  modes  have  opposite 
polarities  (the  sign  of  the  deviation  from  the  constant  background,  see  also  Fig.  7 
below).  The  corresponding  differential  equation  for  $(x)  can  be  solved  in  an  explicit 

form  [cf.  the  above  solution  (9)  for  the  compound  kink)],  yielding  the  following 
expressions  for  the  "larger"  and  "smaller"  bubbles: 


t3  (V j  — t2  )  6  ^*+2/?^^  — i  — (3  (^1  ^2  )  t"sP 


6,  c.  (x,  S)  =  2  tan 

(tj  -t2)  elyx  +  2 p(tx  +t2- 2 t3)eyx  +  p 2 
where  y(S)  =  (l/ 2)  ^(cos£  +  is  +  2)(r3  —  ^ ) (?3 -t2)/ 3  ,  =  -t2) ,  with 


(12) 


r12  3  being  the  roots  of  polynomial  (10),  (1 1 ),  as  they  were  defined  above.  Examples  of 
the  "smaller"  and  "larger"  bubbles  with  opposite  polarities  are  displayed  in  Fig.  7  for 

S  =  0.l7T  . 
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Fig.  7.  (Color  online)  "Smaller"  (a)  and  "larger"  (b)  bubble  solutions,  with  the  opposite 
polarities,  found  at  8  =  O.br . 

Direct  simulations  of  Eqs.  (1)  demonstrate  that  the  bubbles  are  metastable  modes, 
which  persist  for  ~60  time  units,  before  the  numerical  noise  breaks  them  up.  The 
evolution  of  a  smaller  bubble  is  shown  in  Fig.  for  a  relatively  high  damping  parameter, 
a  =  0.5.  The  bubble's  peak  does  not  persist,  relaxing  back  to  to  the  level  corresponding 
to  the  true  vacuum  state.  Thus,  the  bubble  breaks  up  into  kinks  of  types  Ab  and  bA, 
that  very  rapidly  reach  their  equilibrium  velocities,  which  is  considered  in  the  next 
section,  see  Eq.  (21). 
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Fig.  8.  (Color  online)  The  spontaneous  breakup  of  the  "smaller"  bubble  mode  <j>s, 
see  Eq.  (12),  at  a  =  0.5  and  5  =  0.1;r  . 

If  the  dissipation  terms  are  absent,  the  smaller  bubble  breaks  up  into  a  pair  of  Ab 
and  bA  kinks,  that  move  away  from  the  origin  at  a  velocity  approaching  the  Swihart  limit, 
along  with  excitations  in  a  slower  moving  wake,  see  Fig.  9. 


Fig.  9.  (Color  online)  The  breakup  of  the  "smaller"  bubble  mode,  for  the  same  8  =  O.br 
as  in  Fig.  12,  but  without  the  dissipation  (a  =  0 ). 


On  the  other  hand,  if  the  smaller  bubble  evolves  with  initial  non-zero  damping, 
which  is  then  rapidly  reduced  to  zero,  the  evolution  proceeds  differently,  resulting  in  the 
formation  of  a  persistent  breather  at  the  initial  position  of  the  bubble,  see  Fig.  10. 
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Fig.  1 0.  (Color  online)  The  evolution  of  the  "small"  bubble  solution  for  8  =  0.  \n ,  with  the 
damping  switched  off  as  a(r)  =  0.4exp(-r2). 

Finally,  at  values  of  the  dissipation  constant  exceeding  a  certain  critical  value, 
*  0.6 ,  the  bubble's  peak  in  the  corresponding  overdamped  system  collapses  back  to 

the  spatially  uniform  false-vacuum  state  (shallower  potential  minimum),  inhibiting  the 
dynamics  that  would  lead  to  the  emergence  of  the  lowest  potential  minimum,  i.e.,  states 
in  which  the  Ab  and  bA  kinks  expand,  replacing  false  vacuum  with  the  true  one  (not 
shown  here  in  detail). 

The  "larger"  bubble  mode  demonstrates  a  behavior  very  similar  to  that  of  its  "smaller" 
counterpart,  with  the  difference  that  its  rest  mass  is  larger,  allowing  for  more  energy  to 
go  into  the  wake  left  after  the  breakup  of  the  original  bubble.  When  the  damping 
coefficient  is  set  to  zero,  the  larger  bubble  also  breaks  up  into  a  pair  of  Ab  and  bA  kinks, 
that  move  away  from  the  origin  approaching  the  Swihart  velocity,  as  shown  in  Fig.  1 1 . 
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Fig.  1 1 .  (Color  online).  The  evolution  of  the  "larger"  bubble  mode  at  a  =  0  and  8  =  O.br . 

There  is  another  decay  scenario  for  the  larger  bubble,  which  is  not  initiated  by  the 
numerical  noise  alone,  but  may  be  induced  by  a  very  small  relative  displacement  (~10'4) 
of  the  (j>x  and  components.  This  is  shown  in  Fig.  12,  and  additionally  by  means  of  the 

<j\  and  (j>2  profiles  at  r  =  80  in  Fig.  13.  In  terms  of  Fig.  1,  the  final  state  is  a  path  that 

goes  from  b  to  G  along  a  straight  line,  then  to  C  (with  a  bowing  out  of  the  path  toward  B 
as  it  does  so),  and  then  back  to  b  along  a  straight  line. 
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Fig.  12.  (Color  online)  The  evolution  of  components  <j)x  (a)  and  (b),  corresponding  to 
the  second  breakup  scenario  of  the  "larger"  bubble,  for  a  =  0  and  8  =  O.br . 


Fig.  13.  Profiles  of  components  (f>x  and  (the  top  and  bottom  lines,  respectively) 

emerging  from  the  second  breakup  scenario  for  the  "larger"  bubble,  at  a  =  0 
and  8  =  O.br . 


4.  Moving  kinks  (shock  waves)  connecting  true  and  false  vacuum  states 


Kinks  that  connect  true  and  false  vacuum  states,  such  as  corresponding  to 
combinations  fA,  Ab  and  bC  in  Fig.  1,  interpolate  between  local  potential  minima  of 
different  depths,  hence  they  will  be  accelerated  by  the  corresponding  tilted  potential 
[2,6].  In  the  presence  of  the  dissipation,  which  is  accounted  for  coefficient  a  in  Eqs. 
(1),  such  kinks  will  be  accelerated  by  the  tilt  until  they  reach  an  equilibrium  velocity  at 
which  the  tilt-induced  pull  is  balanced  by  the  friction  force.  These  kinks,  which  may  exist 
only  in  the  moving  form,  may  also  be  naturally  considered  as  shock  waves  of  the  phase. 

The  equilibrium  velocity  can  be  predicted  in  an  approximate  analytical  form  by  means 
of  the  momentum-balance  analysis,  similar  to  that  performed  by  McLaughlin  and  Scott 
for  kinks  in  the  single  dc-driven  damped  sine-Gordon  equation  [1,6].  To  this  end,  the 
solution  for  the  kink  of  the  bA  type  (in  terms  of  Fig.  1)  may  be  taken  as  the  solution 
available  at  8  =  0 ,  with  a  shift  to  account  for  the  difference  in  the  boundary  conditions: 

</>f  (x,r)  =  (j)x=0  (x,r)  +  £  /  3, 

where  the  exact  kink  solution  at  zero  frustration  is 


(13) 


Nonlinear  Dynamics  of  Globally  Coupled  Sine  Gordon  Equations 


<j)f  0  (x,r)  =  -2  tan  1 


tanh|(x-vt)/  ^8(l-v2)  j 

s  1 


(14) 


Here  v  is  the  velocity  (v  =  l  corresponds  to  the  Swihart  velocity),  and 

$  {x,t)  =  -<f>*=0{x,t).  (15) 

Note  that  the  substitution  of  Eqs.  (14)  and  (15)  into  expression  (3)  for  the  total 
momentum  (with  r/  =  1),  and  subsequent  expansion  up  to  the  linear  terms  in  the  limit  of 
v  — » 0  makes  it  possible  to  identify  the  effective  mass  of  this  kink  in  the  case  of  <5  =  0 
[5]: 


JbA)  =  6V2(l-V3^/9)«3.355.  (16) 

The  balance  condition  for  momentum  P  implies  the  equailibrium  between  the  driving 
force,  induced  by  the  potential  tilt,  and  friction  force,  induced  by  the  dissipative  term: 

^-  =  -aP^+F«,  =  0-  (17) 

at 

Recall  that  the  general  expression  for  the  total  momentum  of  system  (1 )  is  given  by  Eq. 
(3).  The  consideration  of  derivative  dP /  dt,  in  which  the  time  differentiation  is  performed 
inside  the  integral  expression  (3),  and  the  second  time  derivatives  of  fields  </>l2(x,t)  are 

substituted  according  to  Eqs.  (1),  makes  it  possible  to  identify  the  term  corresponding  to 
the  tilt-induced  driving  force.  Taking  into  regard  the  boundary  conditions  corresponding 
to  the  kink  of  bA  type,  after  somewhat  lengthy  but  straightforward  transformations  the 
driving  force  is  cast  into  the  following  form: 


FlA)  =-3V3sin(£/3).  (18) 

This  force  is  balanced  by  the  frictional  force  in  Eq.  (17)  at  the  equilibrium  value  of  the 
momentum, 


P^  =  -3\f3a  1  sin(c>/3).  (19) 

Finally,  making  use  of  the  relativistic  expression  for  the  kink's  momentum, 

M  w|M)v 

eq  “VT7 


(20) 
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where  the  effective  mass  may  be  taken  per  Eq.  (16),  we  obtain  the  final  result  for  the 
equilibrium  velocity  of  the  moving  kink: 


v 


| bA)  _ 


3V3  sin(^> /3) 

+  27  sin2  (8~/  3) 


(21) 


As  a  typical  example  of  the  corresponding  numerical  simulations  of  Eqs.  (1),  we  can 
take  the  one  for  a  moving  bA  kink  at  8  =  a=  0.\7i .  The  kink,  launched  at  t  =  0  in  the 
form  of  the  explicit  solution  (14),  (15)  for  8  =  a  =  0,  attains  the  equilibrium  velocity 

v^'  =-0.4581,  while  the  analytical  result  for  the  same  parameters  given  by  Eq.  (21)  is 

v^'  =-0.4579.  The  simulations  also  demonstrate  that  the  kink  (shock  wave)  moving  at 
the  equilibrium  velocity  is  completely  stable. 

The  comparison  of  the  shape  of  the  numerical  solution  with  its  analytical 
counterpart  (1 3)-(1 5),  shown  in  Fig.  14,  demonstrates  a  virtually  ideal  agreement,  as 
well  as  for  the  equilibrium  velocity. 


Fig.  14.  (Color  online).  The  analytical  approximation  for  moving  kink  (1 3)-(1 5)  (red)  and 
its  numerically  found  counterpart  (blue),  for  8  =  a  =  0.1;r . 

5.  Collisions  between  moving  kinks 

Collisions  can  occur  between  any  two  moving  kinks  that  share  a  common  true-  or 
false-vacuum  background.  An  example  of  the  former  type  would  be  a  collision  between 
bA  and  Af  kinks  (in  terms  of  Fig.  1),  with  the  Ab-bC  collision  being  an  example  of  the 
latter.  The  tilt  in  the  corresponding  potential  will  act  to  propel  the  kinks  so  as  to 
decrease  the  span  of  the  false  vacuum,  simultaneously  increasing  the  domain  filled  with 
the  true  vacuum. 
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This  also  implies  that  the  kinks  in  the  bA-Ab  configuration  would  move  toward  the 
boundaries,  increasing  the  span  of  true  vacuum  between  them,  while  the  kinks  in  the 
Ab-bA  state  would  move  on  the  collision  course,  to  decrease  the  span  of  false  vacuum 
between  them.  In  the  latter  case,  one  may  expect  that  the  collision  between  the  kinks 
will  lead  to  their  annihilation  into  a  bound  excited  state,  in  the  form  of  a  (slowly 
decaying)  breather.  To  induce  the  kinks  to  move  toward  one  another  In  the  former 
configuration,  it  is  necessary  to  boost  them  in  directions  that  are  opposite  to  those  of 
their  equilibrium  velocities.  In  that  case,  the  kinks  may  rebound,  simultaneously 
generating  radiation,  thus  making  the  collision  complex  (if  the  damping  is  low). 
Examples  of  these  two  types  of  the  collisions  are  presented  below. 

5.1.  Collisions  of  kinks  with  false  vacuums  at  domain  boundaries 


We  first  look  at  the  simplest  collision  of  this  type,  which  is  between  the  kinks  of  the  bA 
and  Ab  types.  There  are  some  practical  constraints  in  carrying  out  the  corresponding 
numerical  experiments,  the  first  of  which  is  that  the  two  kinks  must  be  well  separated 
initially.  Further,  setting  the  two  kinks  far  from  each  other  initially  means  that  they  must 
be  launched  with  sufficiently  large  velocities,  so  that  they  may  collide  before  the  tilt  in 
the  potential  will  reverse  their  trajectories,  as  mentioned  above. 

Another  constraint  for  numerical  simulations  stems  from  the  boundary  conditions  at 
edges  of  the  integration  domain,  x  =  ±L.  Boundary  conditions  of  two  types  have  been 
employed  in  simulations  of  the  collisions:  the  no-radiation  constraint, 


[dtjfjdx  +  d$/Sr]=  0,  x=  -L, 
[0$/cbc  -  0,  x=  +L, 

and  the  zero-magnetic-field  condition, 

[<3^./3x  ]=  0  at  x=  ±L. 


(22) 


(23) 


Since  the  no-radiation  condition  can  still  allow  reflections  of  slowly  moving  disturbances, 
we  have  used,  in  most  cases,  a  combination  of  the  zero-magnetic-field  condition  with  an 
absorber  (a  stripe  with  an  artificially  high  damping  coefficient,  attached  to  each 
boundary). 

The  example  of  the  collision  shown  below  is  for  the  bA-Ab  configuration  with  a  small 
damping  coefficient,  a  =  0.02 ,  and  the  same  value  of  the  frustration  as  used  above,  viz.. 
<5  =  0.1 7:  The  initial  field  was  taken  as  the  joined  pair  of  exact  kinks  for  a  =  5  =  0 ,  with 
the  separation  between  them  2x0  =  30  and  velocities  v  =  ±0.99 . 
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</>l(x,t  =  0)  = 


2  tan 


-1 


tanh^(x-x0)/  j 


a/3(1-v2) 


s 

+  — 
3 


0(x) 


2  tan 


-i 


tanh 


V'x) 

A1-'-2) 


+  - 


0(-*), 


(24) 


as  shown  in  Fig.  15.  a  typical  result  of  the  collision  is  displayed  in  Fig.  16. 


Fig.  15.  (Color  online)  The  initial  bA-Ab  pair  of  boosted  kinks,  prepared  for  the  collision 
the  (j)x  (blue,  with  the  minimum  at  x  =  0 )  and  (red,  with  a  maximum  at  x  =  0 ) 
components. 
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the  collision  of  the  bA-Ab  type,  with  initial  inward  velocities,  v  =  ±0.99  and  damping 
constant  a  =  0.02  . 


In  the  case  of  the  bA-Af  collision,  the  picture  is  essentially  the  same  as  in  Fig.  16,  with 
an  obvious  difference  that  component  <j>2(x,t)  in  the  Af  kink  makes  an  extra  2n  jump. 
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When  the  damping  coefficient  is  set  to  zero,  a  =  0 ,  Fig.  17  shows  that  a  long-lived 
breather  can  be  excited  at  the  origin  by  the  annihilation  of  the  colliding  bA-Af  pair  (which 
is  not  the  case  for  the  bA-Ab  collision).  This  outcome  is  possible  because  the  difference 
between  the  b  and  f  states  in  Fig.  1  is  A</>2  =  2n . 


Fig.  17.  (Color  online)  The  plot  of  </\(x,t)  illustrating  the  formation  of  the  breather 


following  the  bA-Af  collision  with  initial  inward  velocities  v  =  ±0.99  and  zero  damping, 

a  =  0 . 


5.2  .  Collisions  between  traveling  kinks  with  true  vacua  at  boundaries 

Finally,  we  consider  the  collision  of  the  second  type  (as  defined  above),  when  the 
two  kinks  are  supported  by  true  vacuum  states  at  the  spatial  boundaries,  sharing  a  false 
vacuum  state  in  the  region  between  them.  An  example  is  the  collision  between  Ab  and 
bC  kinks,  that  are  launched  toward  the  origin  at  t  =  0  and  at  the  equilibrium  velocities, 
see  Eq.  (21),  from  well-separated  locations,  as  shown  below  in  Fig.  18. 
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(b) 

Fig  18.  (Color  online)  (a)  and  (b):  Plots  of  and  </>2(x,t)  illustrating  in  the  collision 

of  the  Ab  and  bC  kinks  for  a  =  8  =  O.bz- ,  initially  moving  at  the  corresponding  equilibrium 
velocities  [see  Eq.  (21)],  v  =  ±0.4581. 
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As  seen,  this  collision  leads  (quite  naturally)  to  the  fusion  of  the  colliding  kinks  into  a 
quiescent  one  of  the  AC  type.  For  other  pairs  of  moving  kinks  with  the  shared  false 
vacuum  initially  located  between  them,  the  outcome  is  similar,  also  featuring  the  fusion 
into  an  immobile  kink  interpolating  between  two  true  vacuums. 


6.  Conclusions 

This  paper  reports  results  of  the  systematic  analysis  of  the  dynamics  of  fluxons 
(kinks)  in  the  fundamental  prismatic  configuration  formed  by  three  bulk 
superconductors,  which  creates  a  set  of  three  parallel  long  Josephson  junctions  at 
interfaces  between  them.  The  setting  is  described  by  a  system  of  three  coupled  sine- 
Gordon  equations  for  the  phases  corresponding  to  each  junction.  In  fact,  the  condition 
for  the  external  magnetic  flux  trapped  by  the  system  reduces  the  model  to  two 
equations.  The  equations  include  the  dissipative  terms,  and  are  controlled  by  the 
frustration  parameter,  that  measures  the  deviation  of  the  external  flux  from  half  a 
quantum.  Analyzing  the  corresponding  potential  profile,  we  have  identified  different 
types  of  topological  kink  solitons,  which  correspond  to  fluxons,  in  terms  of  the  coupled 
junctions.  Nontopological  "bubble"  modes  have  been  found  too.  Some  solutions, 
including  those  for  "compound"  kinks  and  for  two  types  of  the  bubbles,  were  obtained  in 
an  analytical  form.  Numerical  simulations  demonstrate  that  the  compound  kinks  are 
unstable  against  breaking  up  into  pairs  of  simple  kinks.  Bubbles  are  metastable  objects, 
that  eventually  break  up  into  kink-antikink  pairs.  The  system  also  gives  rise  to  kinks 
which  connect  different  potential  minima,  hence  they  are  pulled  by  the  tilt  of  the 
potential.  Using  the  momentum-balance  method,  we  have  derived  the  equilibrium 
velocity  at  which  such  driven  kinks  should  move,  and  verified  the  prediction  by 
simulations.  Finally,  collisions  between  the  moving  kinks  were  studied  by  means  of 
direct  simulations,  which  demonstrate  that  the  collisions  are  always  strongly  inelastic. 

References 

1.  D.  W.  McLaughlin  and  A.  C.  Scott,  Phys.  Rev.  A  18,  1652  (1978). 

2.  B.  A.  Malomed,  Phys.  Lett.  A  136,  395  (1989). 

3.  G.  S.  Mkrtchyan  and  V.  V.  Shmidt,  Solid  State  Commun.  30,  791  (1979);  B. 
Dueholm,  O.  A.  Levring,  J.  Mygind,  N.  F.  Pedersen,  O.  H.  Soerensen,  and  M. 
Cirillo,  Phys.  Rev.  Lett.  46,  1299  (1981);  E.  Joergensen,  V.  P.  Koshelets,  R. 
Monaco,  J.  Mygind,  M.  R.  Samuelsen,  and  M.  Salerno,  ibid.  49,  1093  (1982);  M. 
Salerno,  M.  R.  Samuelsen,  G.  Filatrella,  S.  Pagano,  and  R.  D.  Parmentier,  Phys. 
Rev.  B  41,  6641  (1990);  J.  G.  Caputo,  N.  Flytzanis,  and  M.  Devorte,  ibid.  50, 
6471  (1994);  R.  G.  Mints  and  I.  B.  Snapiro,  ibid.  52,  9691  (1995);  M.  Cirillo,  N. 
Gronbech-Jensen,  M.  R.  Samuelsen,  M.  Salerno,  and  G.  V.  Rinati,  ibid.  58, 
12377  (1998). 

4.  Yu.  S.  Kivshar  and  B.  A.  Malomed,  Phys.  Rev.  B  37,  9325  (1988);  A.  V.  Ustinov, 
H.  Kohlstedt,  M.  Cirillo,  N.  F.  Pedersen,  G.  Hallmanns,  and  C.  Heiden,  ibid.  B  48, 
10614  (1993);  N.  Gronbech-Jensen  and  M.  R.  Samuelsen,  ibid.  48,  16160 
(1993);  and  J.  A.  Blackburn,  Phys.  Rev.  Lett.  70,  1251  (1993);  N.  Gronbech- 


Nonlinear  Dynamics  of  Globally  Coupled  Sine  Gordon  Equations 


Jensen,  D.  Cai,  A.  R.  Bishop,  A.  W.  G.  Lau,  and  P.  S.  Lomdahl,  Phys.  Rev.  B  50, 
6352  (1994);  Z.  Hermon,  E.  Ben-Jacob,  and  G.  Schon,  ibid.  B  54,  1234  (1996);  S. 
Watanabe,  H.  S.  J.  van  derZant,  S.  H.  Strogatz,  and  T.  P.  Orlando,  Physica  D 
97,  429  (1996);  G.  Hechtfischer,  R.  Kleiner,  K.  Schlenga,  W.  Walkenhorst,  P. 
Muller,  and  H.  L.  Johnson,  Phys.  Rev.  B  55,  14638  (1997);  E.  Goldobin,  A. 

Wall  raff,  N.  Thyssen,  and  A.  V.  Ustinov,  ibid.  57,  130  (1998);  V.  M.  Krasnov,  V.  A. 
Oboznov,  V.  V.  Ryazanov,  N.  Mros,  A.  Yurgens,  and  D.  Winkler,  ibid.  61,  766 
(2000);  E.  Goldobin,  B.  A.  Malomed,  and  A.  V.  Ustinov,  ibid.  62,  1414  (2000);  R. 
Kleiner,  T.  Gaber,  and  G.  Hechtfischer,  ibid.  62,  4086  (2000);  M.  Machida,  T. 
Koyama,  A.  Tanaka,  and  M.  Tachiki,  Physica  C  330,  85  (2000);  E.  Goldobin,  A. 
Wallraff,  and  A.  V.  Ustinov,  J.  Low  Temp.  Phys.  119,  589  (2000);  C.  Gorria,  P.  L. 
Christiansen,  Y.  B.  Gaididei,  V.  Muto,  N.  F.  Pedersen,  and  M.  P.  Soerensen, 
Phys.  Rev.  B  68,  035415  (2003). 

5.  S.  P.  Yukon  and  N.  H.  Lin,  IEEE  Trans.  Mag.  27,  2736  (1991 ). 

6.  Yu.  S.  Kivsharand  B.  A.  Malomed,  Rev.  Mod.  Phys.  61,  763  (1989). 

7.  N.  F.  Pedersen,  Physica  D  68,  27  (1993);  N.  F.  Pedersen  and  A.  V.  Ustinov, 
Supercond.  Science  and  Technology  8,  389  (1995);  A.  V.  Ustinov,  Pysica  D  123, 
315  (1998). 

8.  L.  M.  Floria  and  J.  J.  Mazo,  Adv.  Phys.  45,  505  (1996);  O.  M.  Braun  and  Yu.  S. 
Kivshar,  Phys.  Rep.  306,  1  (1998). 

9.  S.  Savel'ev,  V.  A.  Yampol'skii,  A.  L.  Rakhmanov,  and  F.  Nori,  Rep.  Progr.  Phys. 
73,  026501  (2010). 

10.  J.  Yang  and  S.  P.  Yukon,  unpublished  [AFRL  Final  Report  (F33601-02-F-A581- 
P00010)]. 

11.  T.  Kapitula,  P.  G.  Kevrekidis,  and  Z.  G.  Chen,  SIAM  J.  Appl.  Dyn.  Systems  5,  598 
(2006);  K.  Huybrechts,  G.  Morthier,  and  B.  Maes,  J.  Opt.  Soc.  Am.  B  27,  708 
(2010). 

12.  V.  S.  Shchesnovich  and  D.  S.  Mogilevtsev,  Phys.  Rev.  A  82,  043621  (2010). 

13.  B.  Liu,  L.-B.  Fu,  S.-P.  Yang,  and  J.  Liu,  Phys.  Rev.  A  75,  033601  (2007);  T.  F. 
Viscondi,  K.  Furuya,  and  M.  C.  de  Oliveira,  EPL  90,  10014  (2010). 

14.  T.  Lahaye,  T.  Pfau,  and  L.  Santos,  Phys.  Rev.  Lett.  104,  170404  (2010). 

15.  A.  Gubeskys  and  B.  A.  Malomed,  Eur.  Phys.  J.  D  28,  283  (2004). 

16.  Z.  Yan,  K.  W.  Chow,  and  B.  A.  Malomed,  Chaos,  Solitons  &  Fractals  42,  3013 
(2009). 

17.  A.  Sigler,  B.  A.  Malomed,  and  D.V.  Skryabin,  Phys.  Rev.  E  74,  066604  (2006). 

18.  N.-H.  Kuo,  S.  Sarkar,  and  C.  D.  Hu,  arXiv:0909.1872  (2009). 

19.  I.  V.  Barashenkov  and  E.  Y.  Panova,  Physica  D  69,  114  (1993);  I.  V. 
Barashenkov,  Phys.  Rev.  Lett.  77,  1193  (1996). 


